package com.asa.numanaly;
/**
 * F里面是计算方程组的解，这里也是
 * 只不过这里用迭代法
 * 我对这个不是太感兴趣，所以后两个方法日后实现（主要是我对下一章的问题更感兴趣。。。）
 * p393
 * @author Administrator
 *
 */
public class G {

	
	
	/**
	 * Jacobi迭代法
	 * p393
	 * 这个方法里面的参数，应该还需要更多的自定义
	 * 这个在工程里面是求近似值得，这个特性你懂的，需要定义参数多
	 * 特别的，它这里需要给出近似的x值，而我没有要传，默认定义为都是0了
	 * @param kmax迭代最大次数
	 * @param a
	 * @param b
	 * @return
	 */
	public static double[] fangchengzu1(double a[][],double b[],int kmax) {
		
		
		double q = 0.0001;
		double e = 0.05;
		int i,j,k,n;
		double diag,sum;
		
		
		
		double x[] = {0,0,0,0};
		n = a.length;
		double y[] = new double[n];
		
		for (k = 0; k < kmax; k++) {
			for (int l = 0; l < n; l++) {
				y[l] = x[l];
			}
			for (i = 0; i < n; i++) {
				sum = b[i];
				diag = a[i][i];
				if (Math.abs(diag)<q) {
					System.out.println("diagonal element too small");
					return null;
				}
				for (j = 0; j < n; j++) {
					if (j!=i) {
						sum = sum - a[i][j]*y[j];
					}
				}
				x[i] = sum/diag;
			}
			
			
			for (int l = 0; l < n; l++) {
				System.out.print(x[l]+"\t");
			}
			System.out.println("k:"+k);
			
		}
		double ee = 0;
		for (int l = 0; l < n; l++) {
			ee += Math.sqrt(Math.abs((y[l]-x[l])));
		}
		System.out.println("===========================");
		if (true) {
			for (int l = 0; l < n; l++) {
				System.out.println("k:"+k+"   "+"x"+l+":"+x[l]);
			}
			return null;
			
		}
		
		
		System.out.println("maximum iterations reached");


		
		return x ;


		
	}
	
	
	/**
	 * Gauss-Seidel迭代法
	 * p395		(347,这是另一本书给的，更简洁)
	 * 这个是对Jacob的优化
	 * @param kmax迭代最大次数
	 * @param a
	 * @param b
	 * @return
	 */
	public static double[] fangchengzu2(double a[][],double b[],int kmax) {
		
		double[] x = new double[a.length];
		for (int k = 1; k < kmax; k++) {
			
			for (int i = 0; i < b.length; i++) {
				
				double asa = 0;
				
				for (int j = 0; j < x.length; j++) {
					if (j!=i) {
						asa += a[i][j]*x[j];
					}
				}
				
				x[i] = (b[i]-asa)/a[i][i];
				
				
				
				
			}
			
	
			
		}

		
		return x;
		

		
		
	}
	
	
	
	
	
	/**
	 * SOR算法
	 * p402		
	 * 这个是对Jacob的优化
	 * @param kmax迭代最大次数
	 * @param a
	 * @param b
	 * @param w 0 to 2
	 * @return
	 */
	public static double[] fangchengzu3(double a[][],double b[],int kmax,double w ) {
		
		double[] x = new double[a.length];
		for (int k = 1; k < kmax; k++) {
			
			for (int i = 0; i < b.length; i++) {
				
				double asa = 0;
				
				for (int j = 0; j < x.length; j++) {
					if (j!=i) {
						asa += a[i][j]*x[j];
					}
				}
				
				x[i] = (1-w)*x[i]+ w*(b[i]-asa)/a[i][i];
				//这个地方可以添加精度判断，当前的x与输入的初始值进行比较，用欧氏距离来技术
				
				
				
			}
			
	
			
		}

		
		return x;
		

		
		
	}
	
	
	
	
	/**
	 * 迭代改进
	 * p410	
	 * @param kmax迭代最大次数
	 * @param a
	 * @param b
	 * @param w 0 to 2
	 * @return
	 */
	public static double[] fangchengzu4(double a[][],double b[],int kmax,double w ) {
		
		
		
		
		
		
		return b;
		
		
	}
	
	
	
	
	
	/**
	 * 预处理共轭梯度法
	 * p423
	 * 这个是对Jacob的优化
	 * @param kmax迭代最大次数
	 * @param a
	 * @param b
	 * @param w 0 to 2
	 * @return
	 */
	public static double[] fangchengzu5(double a[][],double b[],int kmax,double w ) {
		
		
		
		
		
		
		return b;
		
		
	}
	
	
	
	
	
	
	
	public static void main(String[] args) {
		double a[][] = {{2,-1,0,0},{-1,2,-1,0},{0,-1,2,-1},{0,0,-1,2}};

		double b[] = {1,0,0,1};

		
		
		
//		fangchengzu1(a, b,100);
		double[] x = fangchengzu2(a, b,100);
		
		for (int i = 0; i < x.length; i++) {
			System.out.println(x[i]);
		}
	}
	
	
	
	
	
	
	
	
	
	
	
}
